docs/rxd-tutorials/Changing initial conditions and parameters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Varying initial concentrations and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A version of this notebook may be run online via Google Colab at https://tinyurl.com/rxd-initial-and-parameters\n",
    " (make a copy or open in playground mode)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often we will want to see how the choice of initial conditions affects the dynamics. We can do this by setting the initial attribute of an `rxd.Species` and rerunning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, suppose at a single point we have the bistable dynamics introduced in the first part of this tutorial. That is, $u'=-u(1-u)(\\alpha - u)$. (Here we use $u$ instead of a specific molecule name to indicate that we are not describing any particular molecule's kinetics, but rather modeling a class of phenomena.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, we'll use an `rxd.Parameter` for $\\alpha$ instead of a constant. This offers two advantages: (1) this allows $\\alpha$ to vary spatially, and (2) this allows us to change the parameter values and rerun without changing the reaction/rate specification."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by defining the model and setting up the recordings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h, rxd\n",
    "from neuron.units import mV, ms, mM\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "h.load_file(\"stdrun.hoc\")\n",
    "\n",
    "soma = h.Section(name=\"soma\")\n",
    "cyt = rxd.Region([soma], name=\"cyt\", nrn_region=\"i\")\n",
    "u = rxd.Species(cyt, name=\"u\")\n",
    "α = rxd.Parameter(cyt, value=0.3 * mM)\n",
    "\n",
    "rate = rxd.Rate(u, -u * (1 * mM - u) * (α - u))\n",
    "\n",
    "h.finitialize(-65 * mV)\n",
    "\n",
    "t = h.Vector().record(h._ref_t)\n",
    "y = h.Vector().record(soma(0.5)._ref_ui)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define a function that sets the initial value for ca, runs the simulation, and plots the curves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_it(ca_init):\n",
    "    u.initial = ca_init\n",
    "    h.finitialize(-65 * mV)\n",
    "    h.continuerun(10 * ms)\n",
    "    plt.plot(t, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now let's run it a few times and plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for u0 in range(-4 * mM, 15 * mM, 2 * mM):\n",
    "    plot_it(u0 * 0.1)\n",
    "\n",
    "plt.xlabel(\"t (ms)\")\n",
    "plt.ylabel(\"[u] (mM)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see the stable fixed points at $0$ and $1$ and how all other trajectories (except for the one starting exactly at $\\alpha$, not shown) tend towards one of those values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we might want to know how the kinetics change if we start at a fixed value (say 0.5) and vary $\\alpha$ instead. We can do this in the same way as above by setting the `value` property of $\\alpha$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_α(α_init):\n",
    "    u.initial = 0.5 * mM\n",
    "    α.value = α_init\n",
    "    h.finitialize(-65 * mV)\n",
    "    h.continuerun(10 * ms)\n",
    "    plt.gca().plot(t, y, label=f\"α = {α_init} mM\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have added a `label` argument to the plot and explicitly specify the `gca()` axes to allow displaying a legend; see the <a href=\"https://matplotlib.org/api/_as_gen/matplotlib.pyplot.legend.html\">matplotlib.pyplot.legend</a> documentation for more.\n",
    "\n",
    "Now we run it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for my_α in [0 * mM, 0.2 * mM, 0.4 * mM, 0.6 * mM, 0.8 * mM, 1 * mM]:\n",
    "    plot_α(my_α)\n",
    "\n",
    "plt.gca().legend()\n",
    "plt.xlabel(\"t (ms)\")\n",
    "plt.ylabel(\"[u] (mM)\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}